library(RColorBrewer); #brewer.pal(n=8, name = "Dark2") #;display.brewer.pal(n = 8, name = "Dark2")
library(tidyverse); theme_set(theme_classic()) #ggplot2, dplyr, tibble, etc.
library(plotly) #use ggplotly() for interactive plots with scoll-over IDs
library(knitr) #use kable() to make formatted tables
library(kableExtra)
library(Rtsne)
library(glmmTMB)
library(lme4)
library(sjPlot)
library(ggpubr) #ggarrange() figure wraps
library(plyr)
library(ggeffects) #ggpredict
library(rstanarm) #stan models
library(shinystan) #stan model evaluation >launch_shinystan_demo()
library(loo) #loo() to compare fits between bayesian models
library(MuMIn) #dredge, model averaging
library(AICcmodavg) #aictab() for AICc model comparisons
library(coin) #permutation test for ratios with independence_test()
library(janitor) #get_dupes() finds and prints duplicates
allAges.mc=tibble(read.csv("CB-mouthColor.csv", h=T)) #N=215
#subset with ages 23-27
mc <- allAges.mc %>% filter(between(age,23,27)) #N=133
Create subsets of mc for females (N=78) and males (N=55) only
#females only
female.mc <- mc %>%
filter(sex == "F") #N=78
#males only
male.mc <- mc %>%
filter(sex == "M") #N=55
Nestling morphometrics are being used here to boost sample size for calculating mass/tarsus residuals
Replace zeros with NAs, sort and filter by age (23-27)
nm=tibble(read.csv("nestlingMorph.csv", h=T))
#replace zeros with NAs in morphometrics
nm[nm==0] <- NA
#remove nestlingMeasurementNumber
nm <- nm[,-1]
#sort by year > nestName > id
nm <- nm %>% arrange(year, nestName ,id)
#filter to include only ages 23-27
nm <- nm %>% filter(between(age, 23, 27))
Create subset of nm including only sexed birds (N=503; females=268; males=235)
#subset nm to include only sexed birds and remove NAs from mass & tarsus
sexed.nm <- nm %>% drop_na(sex, tarsus, mass) #N=503
#make sex a factor rather than character
sexed.nm$sex <- as.factor(sexed.nm$sex)
Only includes sexed individuals (N=503)
sexed.nm %>%
ggplot(aes(mass, fill=sex))+
geom_histogram(binwidth = 10)+
scale_fill_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))+
annotate(geom = "text", x=375, y=29, label="N=268")+
annotate(geom = "text", x=390, y=5, label="N=235")
Only includes sexed individuals (N=503)
sexed.nm %>%
ggplot(aes(tarsus, fill=sex))+
geom_histogram(binwidth = 0.5)+
scale_fill_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))+
annotate(geom = "text", x=59, y=28, label="N=268")+
annotate(geom = "text", x=60, y=5, label="N=235")
massTarsus.plot <- sexed.nm %>%
ggplot(aes(x=tarsus,y=mass,color=sex,label=id))+
geom_point(alpha=0.5)+
geom_smooth()+
scale_color_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))
ggplotly(massTarsus.plot)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Remove outliers (SWB-W BURP03, EN WKAY05, GL JSUP05, YASB WCEM01), now N=499
sexed.nm <- sexed.nm[-c(211, 141, 332, 353),] #N=499
massTarsus.plot <- sexed.nm %>%
ggplot(aes(x=tarsus,y=mass,color=sex,label=id))+
geom_point(alpha=0.5)+
geom_smooth()+
scale_color_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))
ggplotly(massTarsus.plot)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#females only
female.nm <- sexed.nm %>%
filter(sex == "F") #N=265
#males only
male.nm <- sexed.nm %>%
filter(sex == "M") #N=234
femaleIndex.lm <- lm(mass~tarsus, female.nm)
summary(femaleIndex.lm)
##
## Call:
## lm(formula = mass ~ tarsus, data = female.nm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -107.065 -17.796 2.143 18.200 95.956
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -373.9956 44.1717 -8.467 1.8e-15 ***
## tarsus 12.6523 0.7617 16.610 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.89 on 263 degrees of freedom
## Multiple R-squared: 0.512, Adjusted R-squared: 0.5101
## F-statistic: 275.9 on 1 and 263 DF, p-value: < 2.2e-16
Check residual distributions
*Q-Q plot looks a little non-normal, probably seeing the effects of reduced sample size by separating sexes.
plot(femaleIndex.lm)
### Create male mass~tarsus residual index
maleIndex.lm <- lm(mass~tarsus, male.nm)
summary(maleIndex.lm)
##
## Call:
## lm(formula = mass ~ tarsus, data = male.nm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -114.779 -24.276 2.291 26.091 114.680
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -334.1798 52.9473 -6.312 1.39e-09 ***
## tarsus 12.0998 0.8885 13.618 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.43 on 232 degrees of freedom
## Multiple R-squared: 0.4442, Adjusted R-squared: 0.4419
## F-statistic: 185.5 on 1 and 232 DF, p-value: < 2.2e-16
Plot residuals
plot(maleIndex.lm)
Create dataframes containing id and residuals
#male residuals pulled from lm() output
male.nm$resids <- maleIndex.lm$residuals
#female residuals pulled from lm() output
female.nm$resids <- femaleIndex.lm$residuals
#dataframes containing male and female ids and residuals for the join with mc
joinMaleResidual <- data.frame(id = male.nm$id,
resids = male.nm$resids)
joinFemaleResidual <- data.frame(id = female.nm$id,
resids = female.nm$resids)
Join residuals with male and female mouth color data by common id
#left join with female mc
female.mc <- left_join(female.mc,joinFemaleResidual, by="id")
#left join with male mc
male.mc <- left_join(male.mc,joinMaleResidual, by="id")
#left join with both sexes (full mc)
bothSexResid <- data.frame(rbind(joinMaleResidual,joinFemaleResidual))
mc <- left_join(mc,bothSexResid, by="id")
Distribution of residuals for both sexes combined
mc %>%
ggplot(aes(x=resids, fill=sex))+
geom_histogram(binwidth = 10)+
scale_fill_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))
mc$ratio <- mc$mass/mc$tarsus
mc$cubeMass <- mc$mass^(1/3)
mc$cubeRatio <- ((mc$mass^(1/3))/mc$tarsus)*100
independence_test(ratio~as.factor(sex), data = mc)
##
## Asymptotic General Independence Test
##
## data: ratio by as.factor(sex) (F, M)
## Z = -2.6339, p-value = 0.008441
## alternative hypothesis: two.sided
Distribution of mass/tarsus ratio for both sexes combined
mc %>%
ggplot(aes(x=ratio, fill=sex))+
geom_histogram(binwidth = 0.4)+
scale_fill_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))
mc %>%
ggplot(aes(x=mass, fill=sex))+
geom_histogram(binwidth = 20)+
scale_fill_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))
mc %>%
ggplot(aes(x=cubeMass, fill=sex))+
geom_histogram(binwidth = 0.1)+
scale_fill_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))
mc %>%
ggplot(aes(x=cubeRatio, fill=sex))+
geom_histogram(binwidth = 0.3)+
scale_fill_manual(values = c("#1B9E77", "#D73027"),
name = "",
labels = c("Female","Male"))
Plot saturation by age
mc %>%
ggplot(aes(x=sat1,y=age)) +
geom_point()+
geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Plot saturation by mass
satMass.plot <- mc %>%
ggplot(aes(x=sat1,y=mass, label = id)) +
geom_point()+
geom_smooth()
ggplotly(satMass.plot)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Plot age by mass
ageMass.plot <- mc %>%
ggplot(aes(x=age, y=mass, label = id)) +
geom_point()+
geom_smooth()+
geom_boxplot()
ggplotly(ageMass.plot)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Continuous x aesthetic -- did you forget aes(group=...)?
sat1Residual.mdl <- lmer(sat1 ~ age + resids + sex +
bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = mc, REML = FALSE)
summary(sat1Residual.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat1 ~ age + resids + sex + bodTemp1 + ambTemp1 + timeOut1 +
## (1 | family)
## Data: mc
##
## AIC BIC logLik deviance df.resid
## 805.4 828.9 -393.7 787.4 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7578 -0.5450 0.1383 0.5430 2.1740
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 108.51 10.417
## Residual 94.31 9.711
## Number of obs: 100, groups: family, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 126.33319 73.96937 1.708
## age -3.81740 1.54508 -2.471
## resids 0.08675 0.04658 1.862
## sexM -2.14906 2.41129 -0.891
## bodTemp1 1.64250 1.86653 0.880
## ambTemp1 1.30666 0.59721 2.188
## timeOut1 -0.01112 0.14501 -0.077
##
## Correlation of Fixed Effects:
## (Intr) age resids sexM bdTmp1 ambTm1
## age -0.292
## resids 0.133 -0.242
## sexM -0.140 -0.069 0.106
## bodTemp1 -0.836 -0.260 0.020 0.154
## ambTemp1 0.093 0.202 -0.160 0.039 -0.345
## timeOut1 -0.027 -0.088 0.043 0.001 0.066 -0.199
Plot residual model
plot_model(sat1Residual.mdl,
title = "Saturation 1 residual model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
sat1Ratio.mdl <- lmer(sat1 ~ age + ratio + sex +
bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = mc, REML = FALSE)
summary(sat1Ratio.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat1 ~ age + ratio + sex + bodTemp1 + ambTemp1 + timeOut1 + (1 |
## family)
## Data: mc
##
## AIC BIC logLik deviance df.resid
## 805.3 828.8 -393.7 787.3 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6577 -0.5549 0.1515 0.5213 2.0561
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 99.29 9.964
## Residual 97.09 9.853
## Number of obs: 100, groups: family, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 93.6305265 73.4209470 1.275
## age -3.6954789 1.4948678 -2.472
## ratio 5.0348153 2.5351722 1.986
## sexM -3.8538809 2.4927588 -1.546
## bodTemp1 1.5537491 1.8767843 0.828
## ambTemp1 1.4185318 0.5787127 2.451
## timeOut1 0.0005294 0.1435943 0.004
##
## Correlation of Fixed Effects:
## (Intr) age ratio sexM bdTmp1 ambTm1
## age -0.220
## ratio -0.100 -0.213
## sexM -0.127 0.006 -0.240
## bodTemp1 -0.851 -0.266 0.008 0.146
## ambTemp1 0.139 0.190 -0.082 0.078 -0.356
## timeOut1 -0.039 -0.100 0.094 -0.031 0.062 -0.206
Plot ratio model
sat1RatioMdl.plot <- plot_model(sat1Ratio.mdl,
title = "Saturation 1 ratio model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
sat1RatioMdl.plot
sat1Mass.mdl <- lmer(sat1 ~ age + mass + sex +
bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = mc, REML = FALSE)
summary(sat1Mass.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat1 ~ age + mass + sex + bodTemp1 + ambTemp1 + timeOut1 + (1 |
## family)
## Data: mc
##
## AIC BIC logLik deviance df.resid
## 806.4 829.9 -394.2 788.4 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6188 -0.5994 0.1448 0.5711 1.9641
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 102.29 10.114
## Residual 97.53 9.876
## Number of obs: 100, groups: family, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 96.38436 73.84681 1.305
## age -3.48163 1.49686 -2.326
## mass 0.05982 0.03622 1.651
## sexM -4.18234 2.59786 -1.610
## bodTemp1 1.55714 1.88466 0.826
## ambTemp1 1.48949 0.58264 2.556
## timeOut1 0.00317 0.14519 0.022
##
## Correlation of Fixed Effects:
## (Intr) age mass sexM bdTmp1 ambTm1
## age -0.233
## mass -0.097 -0.164
## sexM -0.111 0.015 -0.355
## bodTemp1 -0.849 -0.267 0.008 0.140
## ambTemp1 0.130 0.176 -0.018 0.062 -0.354
## timeOut1 -0.042 -0.100 0.120 -0.050 0.063 -0.200
Plot mass model
sat1MassMdl.plot <- plot_model(sat1Mass.mdl,
title = "Saturation 1 mass model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
sat1MassMdl.plot
sat1CubeMass.mdl <- lmer(sat1 ~ age + cubeMass + sex +
bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = mc, REML = FALSE)
summary(sat1CubeMass.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat1 ~ age + cubeMass + sex + bodTemp1 + ambTemp1 + timeOut1 +
## (1 | family)
## Data: mc
##
## AIC BIC logLik deviance df.resid
## 806.2 829.7 -394.1 788.2 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6156 -0.5987 0.1421 0.5741 1.9641
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 100.89 10.044
## Residual 97.72 9.886
## Number of obs: 100, groups: family, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.329008 80.194046 0.653
## age -3.465349 1.487939 -2.329
## cubeMass 9.462185 5.474152 1.729
## sexM -4.259462 2.598636 -1.639
## bodTemp1 1.500293 1.884291 0.796
## ambTemp1 1.503016 0.580330 2.590
## timeOut1 0.005269 0.144882 0.036
##
## Correlation of Fixed Effects:
## (Intr) age cubMss sexM bdTmp1 ambTm1
## age -0.164
## cubeMass -0.403 -0.154
## sexM 0.009 0.012 -0.356
## bodTemp1 -0.778 -0.266 -0.008 0.145
## ambTemp1 0.123 0.176 -0.008 0.059 -0.355
## timeOut1 -0.078 -0.099 0.125 -0.052 0.060 -0.199
Plot cube mass model
sat1CubeMassMdl.plot <- plot_model(sat1CubeMass.mdl,
title = "Saturation 1 cube mass model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
sat1CubeMassMdl.plot
sat1CubeRatio.mdl <- lmer(sat1 ~ age + cubeRatio + sex +
bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = mc, REML = FALSE)
summary(sat1CubeRatio.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat1 ~ age + cubeRatio + sex + bodTemp1 + ambTemp1 + timeOut1 +
## (1 | family)
## Data: mc
##
## AIC BIC logLik deviance df.resid
## 806.2 829.7 -394.1 788.2 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8199 -0.5831 0.1268 0.5632 2.2226
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 120.82 10.992
## Residual 91.71 9.576
## Number of obs: 100, groups: family, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 55.83330 81.10808 0.688
## age -3.70334 1.58647 -2.334
## cubeRatio 5.55821 3.57531 1.555
## sexM -2.00929 2.40637 -0.835
## bodTemp1 1.63945 1.85847 0.882
## ambTemp1 1.27632 0.61714 2.068
## timeOut1 -0.01771 0.14746 -0.120
##
## Correlation of Fixed Effects:
## (Intr) age cubeRt sexM bdTmp1 ambTm1
## age -0.177
## cubeRatio -0.413 -0.210
## sexM -0.202 -0.074 0.152
## bodTemp1 -0.763 -0.248 0.002 0.152
## ambTemp1 0.163 0.194 -0.180 0.024 -0.329
## timeOut1 -0.036 -0.077 0.003 0.002 0.070 -0.186
Plot cube mass model
sat1CubeRatioMdl.plot <- plot_model(sat1CubeRatio.mdl,
title = "Saturation 1 cube ratio model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
sat1CubeRatioMdl.plot
### Now let’s compare the fit of models with different condition indices
indices.mdls <- c(sat1Residual.mdl, sat1Ratio.mdl, sat1Mass.mdl, sat1CubeMass.mdl, sat1CubeRatio.mdl)
indices.modnames <- c("Residual", "Ratio", "Mass", "Cube Mass", "Cube Ratio")
confset(indices.mdls, modnames = indices.modnames, method = "ordinal")
##
## Confidence set for the best model
##
## Method: ordinal ranking based on delta AIC
##
## Models with substantial weight:
## K AICc Delta_AICc AICcWt
## Ratio 9 807.32 0.00 0.26
## Residual 9 807.43 0.11 0.25
## Cube Mass 9 808.22 0.89 0.17
## Cube Ratio 9 808.22 0.90 0.17
## Mass 9 808.42 1.09 0.15
##
##
## Models with some weight:
## K AICc Delta_AICc AICcWt
##
##
## Models with little weight:
## K AICc Delta_AICc AICcWt
##
##
## Models with no weight:
## K AICc Delta_AICc AICcWt
Using AICc Weight (AICcWt), the confset() function calculates the probability of each model given the data and the other candidate models. So in this case, we have more confidence in the ratio model’s ability to explain the variation in the data, even though AICc values are only slightly different.
hue1.mdl <- lmer(hue1*100 ~ age + ratio + sex + bodTemp1 + ambTemp1 + (1|family),
data = mc, REML = FALSE)
hue1Mdl.plot <- plot_model(hue1.mdl,
title = "hue 1 model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
hue1Mdl.plot
#create a time 2 data frame (N=103)
time2.mc <- mc %>% select(row:sex, ends_with("2"), ambFinal, broodSize, date:skull, resids, ratio)
#I'm dropping ambTemp2 from data frame before dropping NAs.
#ambTemp2 has many more NAs. We can use ambFinal instead.
#This boosts sample size from 32 to 99 obs.
time2.mc <- time2.mc %>% select(!ambTemp2)
#Now drop NAs after removing ambTemp2
time2.mc <- drop_na(time2.mc)
At time 1, we have 78 females and 55 males. At time 2, we have 64 females and 39 males.
sat2.mdl <- lmer(sat2 ~ age + ratio + sex + bodTemp2 + ambFinal + timeOut2 + (1|family),
data = time2.mc, REML = FALSE)
summary(sat2.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat2 ~ age + ratio + sex + bodTemp2 + ambFinal + timeOut2 + (1 |
## family)
## Data: time2.mc
##
## AIC BIC logLik deviance df.resid
## 834.5 858.2 -408.2 816.5 94
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.81870 -0.61972 0.01104 0.52422 2.44593
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 144.83 12.034
## Residual 89.75 9.474
## Number of obs: 103, groups: family, 38
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 84.243482 66.498688 1.267
## age -3.949011 1.631450 -2.421
## ratio 1.828675 2.496676 0.732
## sexM 4.049973 2.466390 1.642
## bodTemp2 2.230622 1.679602 1.328
## ambFinal 1.806990 0.567272 3.185
## timeOut2 0.005522 0.064403 0.086
##
## Correlation of Fixed Effects:
## (Intr) age ratio sexM bdTmp2 ambFnl
## age -0.336
## ratio -0.004 -0.124
## sexM -0.021 0.038 -0.378
## bodTemp2 -0.778 -0.258 -0.168 0.064
## ambFinal 0.047 0.049 0.095 0.083 -0.257
## timeOut2 -0.101 -0.037 0.137 -0.129 0.057 -0.097
Plot sat2 model
sat2Mdl.plot <- plot_model(sat2.mdl,
title = "Saturation 2 ratio model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
sat2Mdl.plot
hue2.mdl <- lmer(hue2*100 ~ age + ratio + sex + bodTemp2 + ambFinal + timeOut2 + (1|family),
data = time2.mc, REML = FALSE)
summary(hue2.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: hue2 * 100 ~ age + ratio + sex + bodTemp2 + ambFinal + timeOut2 +
## (1 | family)
## Data: time2.mc
##
## AIC BIC logLik deviance df.resid
## 565.2 588.9 -273.6 547.2 94
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.47645 -0.47349 0.00881 0.53622 2.39858
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 6.795 2.607
## Residual 7.731 2.780
## Number of obs: 103, groups: family, 38
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 114.52826 17.46943 6.556
## age 0.00711 0.38966 0.018
## ratio -1.63567 0.66936 -2.444
## sexM 0.03974 0.69593 0.057
## bodTemp2 -0.87323 0.46688 -1.870
## ambFinal -0.27935 0.13579 -2.057
## timeOut2 -0.04941 0.01737 -2.845
##
## Correlation of Fixed Effects:
## (Intr) age ratio sexM bdTmp2 ambFnl
## age -0.210
## ratio -0.001 -0.144
## sexM -0.022 0.033 -0.338
## bodTemp2 -0.826 -0.306 -0.161 0.055
## ambFinal 0.097 0.080 0.101 0.107 -0.290
## timeOut2 -0.115 -0.056 0.094 -0.115 0.084 -0.124
Plot hue 2 model
hue2Mdl.plot <- plot_model(hue2.mdl,
title = "hue 2 model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
hue2Mdl.plot
Wrap sat and hue model plots
ratioMdlWrap <- ggarrange(sat1RatioMdl.plot, hue1Mdl.plot, sat2Mdl.plot, hue2Mdl.plot)
ratioMdlWrap
#ggsave("ratioMdlWrap.png", plot = ratioMdlWrap, device = "png", dpi = 300, width = 10, height = 8)
The semi_join() function returns all rows from mc where there are matching IDs in time2.mc, keeping just columns from mc.
bothTimes.mc <- semi_join(mc,time2.mc, by="id") #N=103
sat1.BT.mdl <- lmer(sat1 ~ age + ratio + sex +
bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = bothTimes.mc, REML = FALSE)
summary(sat1.BT.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat1 ~ age + ratio + sex + bodTemp1 + ambTemp1 + timeOut1 + (1 |
## family)
## Data: bothTimes.mc
##
## AIC BIC logLik deviance df.resid
## 715.8 738.2 -348.9 697.8 80
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0229 -0.5839 0.1569 0.6019 1.8685
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 122.56 11.071
## Residual 85.07 9.223
## Number of obs: 89, groups: family, 33
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 86.87942 77.80080 1.117
## age -3.81232 1.60012 -2.383
## ratio 3.86820 2.67839 1.444
## sexM -3.62263 2.51583 -1.440
## bodTemp1 2.05087 1.95853 1.047
## ambTemp1 1.26256 0.63269 1.996
## timeOut1 -0.01961 0.15890 -0.123
##
## Correlation of Fixed Effects:
## (Intr) age ratio sexM bdTmp1 ambTm1
## age -0.245
## ratio -0.104 -0.231
## sexM -0.143 0.032 -0.259
## bodTemp1 -0.847 -0.249 0.028 0.152
## ambTemp1 0.106 0.187 -0.101 0.099 -0.330
## timeOut1 0.006 -0.099 0.053 -0.055 0.030 -0.203
Plot both times sat1 model
sat1.BT.mdl.plot <- plot_model(sat1.BT.mdl,
title = "Both times sat1 model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
sat1.BT.mdl.plot
hue1.BT.mdl <- lmer(hue1*100 ~ age + ratio + sex + bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = bothTimes.mc, REML = FALSE)
Plot both times hue1 model
hue1.BT.mdl.plot <- plot_model(hue1.BT.mdl,
title = "Both times hue1 model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
hue1.BT.mdl.plot
Wrap both times model plots
bothTimes.mdlWrap <- ggarrange(sat1.BT.mdl.plot, hue1.BT.mdl.plot, sat2Mdl.plot, hue2Mdl.plot)
bothTimes.mdlWrap
#ggsave("bothTimesMdlWrap.png", plot = bothTimes.mdlWrap, device = "png", dpi = 300, width = 10, height = 8)
Plot saturation by hue
mc %>% ggplot(aes(x=sat1, y=hue1))+
geom_jitter()+
geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Mean yellow 1 model
meanY.mdl <- lmer(meanY1 ~ age + ratio + sex +
bodTemp1 + ambTemp1 + (1|family),
data = mc, REML = FALSE)
Plot mean yellow 1 model
meanY.plot <- plot_model(meanY.mdl,
title = "mean yellow 1 model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
Mean magenta 1 model
meanM.mdl <- lmer(meanM1 ~ age + ratio + sex +
bodTemp1 + ambTemp1 + (1|family),
data = mc, REML = FALSE)
Plot mean magenta 1 model
meanM.plot <- plot_model(meanM.mdl,
title = "mean magenta 1 model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
Wrap plot of mean Y1 & mean M1 models
ggarrange(meanY.plot,meanM.plot)
Mean yellow 2 model
meanY2.mdl <- lmer(meanY2 ~ age + ratio + sex +
bodTemp2 + ambFinal + timeOut2 + (1|family),
data = time2.mc, REML = FALSE)
Plot mean yellow 2 model
meanY2.plot <- plot_model(meanY2.mdl,
title = "mean yellow 2 model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
Mean magenta 2 model
meanM2.mdl <- lmer(meanM2 ~ age + ratio + sex +
bodTemp2 + ambFinal + timeOut2 + (1|family),
data = time2.mc, REML = FALSE)
Plot mean magenta 2 model
meanM2.plot <- plot_model(meanM2.mdl,
title = "mean magenta 2 model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
Wrap plots of mean Y2 & M2 models
ggarrange(meanY2.plot,meanM2.plot)
meanColorMdlWrap <- ggarrange(meanY.plot,meanM.plot,hue1Mdl.plot, meanY2.plot,meanM2.plot,hue2Mdl.plot)
meanColorMdlWrap
#ggsave("meanColorMdlWrap.png", plot = meanColorMdlWrap, device = "png", dpi = 300, width = 10, height = 8)
Data frame including birds 23-33 days old
older.mc <- allAges.mc %>% filter(between(age,23,33))
#create mass/tarsus ratio column
older.mc$ratio <- older.mc$mass/older.mc$tarsus
Mean yellow 1 model with older birds
olderMeanY.mdl <- lmer(meanY1 ~ age + ratio + sex +
bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = older.mc, REML = FALSE)
Plot mean yellow 1 model with older birds
olderMeanY.plot <- plot_model(olderMeanY.mdl,
title = "mean yellow model with older birds",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
Mean magenta 1 model with older birds
olderMeanM.mdl <- lmer(meanM1 ~ age + ratio + sex +
bodTemp1 + ambTemp1 + (1|family),
data = older.mc, REML = FALSE)
Plot mean magenta 1 model with older birds
olderMeanM.plot <- plot_model(olderMeanM.mdl,
title = "mean magenta model with older birds",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
Wrap mean colors with older birds models
ggarrange(olderMeanY.plot,olderMeanM.plot,meanY.plot,meanM.plot)
#define bin breaks
binBreaks <- c(22,25,28,33) #define breaks
#name bins
bins <- c("bin1","bin2","bin3") #label bins
#assign age variable to bins
older.mc$ageBins <- cut(older.mc$age,breaks = binBreaks,labels = bins) #vectorize
#create a new column for age bins
older.mc$ageBins <- as.factor(older.mc$ageBins)
#checkBins <- data.frame(older.mc$id, older.mc$age, older.mc$ageBins)
ageBinsMeanY.mdl <- lmer(meanY1 ~ ageBins + ratio + sex +
bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = older.mc, REML = FALSE)
ageBinsMeanY.plot <- plot_model(ageBinsMeanY.mdl,
title = "age bin mean yellow model",
bpe = "median",
vline.color = "black",
show.values = TRUE,
value.offset = .3,
line.size = 1.2)
ggarrange(olderMeanY.plot, ageBinsMeanY.plot)
Birds in the oldest category (29-33) are driving the negative age effect.
Comparison of model fit for continuous age and age bin models
ageBins.mdls <- c(olderMeanY.mdl, ageBinsMeanY.mdl)
ageBins.modnames <- c("continuous age", "age bins")
confset(ageBins.mdls, modnames = ageBins.modnames, method = "ordinal")
##
## Confidence set for the best model
##
## Method: ordinal ranking based on delta AIC
##
## Models with substantial weight:
## K AICc Delta_AICc AICcWt
## continuous age 9 1021.89 0 0.95
##
##
## Models with some weight:
## K AICc Delta_AICc AICcWt
## age bins 10 1027.84 5.95 0.05
##
##
## Models with little weight:
## K AICc Delta_AICc AICcWt
##
##
## Models with no weight:
## K AICc Delta_AICc AICcWt
There’s no difference in the model probabilities.